clear all; clc;close all;
%Updated: 05/07/2018

a   =load('agrid2.txt');
b   =load('bgrid2.txt');

vcon0=load('gvcon.txt');
vaut0=load('gvaut.txt');

p0  =load('gp.txt');
h0  =load('gh.txt');
ibp0=load('gibp.txt');
ct0 =load('gct.txt');
ca0 =load('gca.txt');
gn0 =load('ggn.txt');
tau0=load('gtau.txt');
w0  =load('gw.txt');
def0=load('gdef1.txt');
q0  =load('gq.txt');

paut0  =load('gpaut.txt');
haut0  =load('ghaut.txt');
ctaut0 =load('gctaut.txt');
gnaut0 =load('ggnaut.txt');
tauaut0=load('gtauaut.txt');
waut0  =load('gwaut.txt');

na=rows(a);
nk=1;
ns=na/nk;
nb=rows(b);

alpha = 0.75;
pm     = 1;
gamma  = 0.43;
omegac = 0.3;
mu     = 1;
lambda = 0.184;
r      = 0.02;

atrans = ones(nb,1)*a';

for i=1:na;
    for j=1:nb;
        vcon(i,j)  =vcon0((i-1)*nb + j);
        
        p(i,j)  =p0((i-1)*nb + j);
        h(i,j)  =h0((i-1)*nb + j);
        ibp(i,j)=ibp0((i-1)*nb + j);
        ct(i,j) =ct0((i-1)*nb + j);
        ca(i,j) =ca0((i-1)*nb + j);
        gn(i,j) =gn0((i-1)*nb + j);       
        tau(i,j)=tau0((i-1)*nb + j);       
        w(i,j)  =w0((i-1)*nb + j);   
        def(i,j)=def0((i-1)*nb + j);   
        q(i,j)  =q0((i-1)*nb + j);   
        rec(i,j)=w(i,j)*h(i,j)*tau(i,j);
        
        ratio(i,j)=gn(i,j)/(h(i,j)^alpha);   

        if (def(i,j)==0);ratio(i,j)=NaN;end;
    end;
end;

bp = b(ibp);

for i=1:na;
    vaut(i)  =vaut0(i);
    
    paut(i)  =paut0(i);
    haut(i)  =haut0(i);
    ctaut(i) =ctaut0(i);
    gnaut(i) =gnaut0(i);       
    tauaut(i)=tauaut0(i);       
    waut(i)  =waut0(i);   
    recaut(i)=waut(i)*haut(i)*tauaut(i);
	wwaut(i) = alpha*paut(i).*haut(i).^(alpha-1); 
    
end;

cn       = h.^alpha-gn;
cnaut    = haut.^alpha-gnaut;
const    = (omegac.*ct.^(-mu)+(1-omegac).*cn.^(-mu)).^(-1/mu);
constaut = (omegac.*ctaut.^(-mu)+(1-omegac).*cnaut.^(-mu)).^(-1/mu);

for i=1:na;
    temp1 = find(ct(i,:) > 0);
    temp2 = find(cn(i,:) > 0);
    temp3 = find(p(i,:) > 0);

    iminct(i) = temp1(1);
    imincn(i) = temp2(1);
    iminp(i) = temp3(1);
    
    imin(i)  = max(iminct(i),max(imincn(i),iminp(i)));
    
    v(i,1:imin(i)-1)  = NaN; 
    p(i,1:imin(i)-1)  = NaN; 
    bp(i,1:imin(i)-1) = NaN; 
    ibp(i,1:imin(i)-1)= NaN; 
    ctE(i,1:imin(i)-1) = NaN; 
    ctU(i,1:imin(i)-1) = NaN; 
    cnE(i,1:imin(i)-1) = NaN; 
    cnU(i,1:imin(i)-1) = NaN; 
    constE(i,1:imin(i)-1) = NaN;     
    constU(i,1:imin(i)-1) = NaN;     
    gn(i,1:imin(i)-1) = NaN; 
    h(i,1:imin(i)-1)  = NaN; 
    ca(i,1:imin(i)-1)  = NaN; 
    tau(i,1:imin(i)-1)= NaN; 
	rec(i,1:imin(i)-1)= NaN;
	ww(i,1:imin(i)-1)= NaN;
	q(i,1:imin(i)-1)  = NaN;
  
end;

for i=1:ns;
for j=1:nk;        
    vcon2(i,j,:)  = vcon((i-1)*nk+j,:); 
    p2(i,j,:)  = p((i-1)*nk+j,:);  
    bp2(i,j,:) = bp((i-1)*nk+j,:);  
    ibp2(i,j,:)= ibp((i-1)*nk+j,:);  
    ct2(i,j,:) = ct((i-1)*nk+j,:);  
    cn2(i,j,:) = cn((i-1)*nk+j,:);  
    const2(i,j,:) = const((i-1)*nk+j,:);         
    gn2(i,j,:) = gn((i-1)*nk+j,:);  
	ca2(i,j,:) = ca((i-1)*nk+j,:);  
    h2(i,j,:)  = h((i-1)*nk+j,:);   
    tau2(i,j,:)= tau((i-1)*nk+j,:);  
	rec2(i,j,:)= rec((i-1)*nk+j,:); 
    def2(i,j,:)= def((i-1)*nk+j,:);  
	q2(i,j,:)  = q((i-1)*nk+j,:); 
	ww2(i,j,:) = alpha*p2(i,j,:).*h2(i,j,:).^(alpha-1); 
    

    vaut2(i,j)  = vaut((i-1)*nk+j); 
    paut2(i,j)  = paut((i-1)*nk+j);  
    ctaut2(i,j) = ctaut((i-1)*nk+j);  
    cnaut2(i,j) = cnaut((i-1)*nk+j);  
    constaut2(i,j) = constaut((i-1)*nk+j);          
    gnaut2(i,j) = gnaut((i-1)*nk+j);  
    haut2(i,j)  = haut((i-1)*nk+j);  
    tauaut2(i,j)= tauaut((i-1)*nk+j);  
	wwaut2(i,j) = alpha*paut2(i,j).*haut2(i,j).^(alpha-1); 
end;
end;

nettax2    = p2.*gn2-ca2;           %net taxes
nettaxaut2 = paut2.*gnaut2;         %net taxes

fiscal2 = tau2 - p2.*gn2;           %fiscal surplus
fiscalaut2=tauaut2-paut2.*gnaut2;   %fiscal surplus

nfiscal2 = nettax2 - p2.*gn2;           %fiscal surplus net of transfers
nfiscalaut2=nettaxaut2-paut2.*gnaut2;   %fiscal surplus net of transfers

lowa = round(2*ns/6);
meda = round(ns/2);
hgha = round(4*ns/6);

b =-b/(lambda+r);  %change of sign

a = exp(a);
%debt indices
%default IR: 
%%=122 mean debt, =196 0.5*mean debt,   =159 0.75*mean debt, 

%flex wages IR:
%%=113 mean debt, =169 0.5*mean debt, =141 0.75*mean debt,  

%risk-free debt IR:
%%=1888 mean debt, =1944 0.5*mean debt, =1916 0.75*mean debt,  =1774 2*meandebt 
%=1973 0.25*mean debt, =2001 zero, 

lowb = 69;     
hghb = 13;      %mean debt

%put together policy fcns of repayment and autarky
mat_ones = ones(size(def));
risk_free_debt=0;

if (risk_free_debt==1);
    p3  =p;
    h3  =h;
    ct3=ct;
    cn3=cn;
    ca3 =ca;
    gn3 =gn;       
    tau3=tau;       
    w3  =w;   
    rec3=rec;
else;
    p3  =def.*p+(mat_ones-def).*(paut'*ones(1,nb));
    h3  =def.*h+(mat_ones-def).*(haut'*ones(1,nb));
    ct3=def.*ct+(mat_ones-def).*(ctaut'*ones(1,nb));
    cn3=def.*cn+(mat_ones-def).*(cnaut'*ones(1,nb));
    ca3 =def.*ca+(mat_ones-def).*0;
    gn3 =def.*gn+(mat_ones-def).*(gnaut'*ones(1,nb));       
    tau3=def.*tau+(mat_ones-def).*(tauaut'*ones(1,nb));       
    w3  =def.*w+(mat_ones-def).*(waut'*ones(1,nb));   
    rec3=def.*rec+(mat_ones-def).*(recaut'*ones(1,nb));
end;

%cnaut3 = haut3.*cnEaut3 + (1-haut3).*cnUaut3;
%ctaut3 = haut3.*ctEaut3 + (1-haut3).*ctUaut3;

[i1,i2] = find(def>0);
q3  =NaN; bp3 =NaN; sp3 =NaN; gn3 = NaN; h3=NaN; p3=NaN;

for i=1:na;
for j=1:nb;
    if (def(i,j)>0); q3(i,j)=q(i,j);bp3(i,j)=bp(i,j);sp3(i,j)=100*((1+q3(i,j).^(-1)-lambda)./(1+r)-1);p3=p;h3=h;gn3=gn;tau3=tau;w3=w;rec3=rec;
    else;
        q3(i,j)=NaN;bp3(i,j)=NaN;sp3(i,j)=NaN;p3=NaN;h3=NaN;gn3=NaN;tau3=NaN;w3=NaN;rec3=NaN;
    end;
end;
end;

nettax3    = p3.*gn3-ca3;           %net taxes
fiscal3 = tau3 - p3.*gn3;           %fiscal surplus
nfiscal3 = nettax3 - p3.*gn3;           %fiscal surplus net of transfers

%default decision and prices
def21(:,:) = def2(:,1,:);
def22(:,:) = def2(:,1,:);

price21(:,:) = q2(:,1,:);
price22(:,:) = q2(:,1,:);


% figure;
% plot(b(1:nb),ctE21(lowa,:),'-r',b(1:nb),ctE21(meda,:),':r',b(1:nb),ctE21(hgha,:),'--r',...
%     b(1:nb),ctU21(lowa,:),'-b',b(1:nb),ctU21(meda,:),':b',b(1:nb),ctU21(hgha,:),'--b','LineWidth',1.5);
% xlabel('b');legend('c^T_E y^T low','c^T_E y^T med','c^T_E y^T hgh','c^T_U y^T low','c^T_U y^T med','c^T_U y^T hgh');
% xlim([min(b) 1]);
% 
% figure;
% subplot(2,2,1);
% contourf(b(1:nb),a(1:ns),def21(:,:),20);
% xlabel('b');
% ylabel('a');
% title('default decision \kappa_{L}');
% subplot(2,2,2);
% contourf(b(1:nb),a(1:ns),def22(:,:),20);
% xlabel('b');
% ylabel('a');
% title('default decision \kappa_{H}');
% subplot(2,2,[3 4]);
% plot(b(1:nb),price21(lowa,:),'-r',b(1:nb),price21(meda,:),':r',b(1:nb),price21(hgha,:),'--r',...
%     b(1:nb),price22(lowa,:),'-b',b(1:nb),price22(meda,:),':b',b(1:nb),price22(hgha,:),'--b','LineWidth',1.5);
% xlabel('b');
% 
% 
% %%
% %value functions: repayment vs autarky
% vcon21(:,:) = vcon2(:,1,:);
% vcon22(:,:) = vcon2(:,1,:);
% vaut21(:,:) = vaut2(:,1);
% vaut22(:,:) = vaut2(:,1);
% 
% figure;
% plot(b(1:nb),vcon21(lowa,:),'-r',b(1:nb),vcon21(meda,:),':r',b(1:nb),...
%     vaut21(lowa,:)*ones(size(b)),'-b',b(1:nb),vaut21(meda)*ones(size(b)),':b','LineWidth',1.5);
% xlabel('b');title('value functions');
% ylim([-16 -13]);legend('rep low y^T','rep high y^T','aut low y^T','aut high y^T')

%% Plot variables as function of yT for two debt levels
%optimal policies and prices: repayment vs autarky
%NameArray = {'LineStyle','Color','LineWidth'};
%ValueArray = {'-',':','-',':';'red','red','blue','blue';1.5,1.5,1.5,1.5}';

%NameArray0 = {'LineStyle','Color','LineWidth'};
%ValueArray0 = {'-','-';'red','blue';1.5,1.5}';
NameArray = {'LineStyle','Color','LineWidth'};
NameArray0 = {'LineStyle','Color','LineWidth'};
ValueArray0 = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';
ValueArray = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';
ValueArray2 = {'--','-','--','-';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray4 = {'--','--','-','-';'red','blue','red','blue';1.5,1.5,1.5,1.5}';
ValueArray5 = {':',':';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';

endb = 2;
lowb = 69;
hghb = 13;

figure;
subplot(5,2,1);P = plot(a,100*squeeze(1-h3(:,lowb)),a,100*squeeze(1-h3(:,hghb)));
title('unemployment','Fontweight','normal');
AX1 = gca;box on;ylabel('%');
set(AX1,'YTick',0:10:20);
axis(AX1, [a(1) a(end) -1 20]);
set(P,NameArray,ValueArray);
subplot(5,2,2);P = plot(a,squeeze(p3(:,lowb)),a,squeeze(p3(:,hghb)));
title('p^{N}','Fontweight','normal'); 
%lgd=legend('low y^T','high y^T');lgd.Location='SouthEast';
legend('low debt (0.5b^{ss})','high debt (b^{ss})','Location','SouthEast','Orientation','vertical');
AX1 = gca;box on;
set(AX1,'YTick',1:2:5);
axis(AX1, [a(1) a(end) 1 5]);
set(P,NameArray,ValueArray);
subplot(5,2,3);P = plot(a,squeeze(gn3(:,lowb)),a,squeeze(gn3(:,hghb)));
title('g^{N}','Fontweight','normal'); 
AX1 = gca;box on;
set(AX1,'YTick',0:0.2:0.4);
axis(AX1, [a(1) a(end) 0 0.4]);
set(P,NameArray,ValueArray);
%subplot(5,2,4);P = plot(b,squeeze(bpr2(lowa,1,:)),b,squeeze(bpr2(hgha,1,:)));
subplot(5,2,4);P = plot(a,-squeeze(bp3(:,lowb))/(lambda+r),a,-squeeze(bp3(:,hghb))/(lambda+r));
title('debt','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.25:0.5:1.25);
axis(AX1, [a(1) a(end) 0.25 1.25]);
set(P,NameArray0,ValueArray0);
subplot(5,2,5);P = plot(a,squeeze(cn3(:,lowb)),a,squeeze(cn3(:,hghb)));
title('c^{N}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.8:0.05:0.9);
axis(AX1, [a(1) a(end) 0.8 0.9]);
set(P,NameArray,ValueArray);
subplot(5,2,6);P = plot(a,squeeze(ct3(:,lowb)),a,squeeze(ct3(:,hghb)));
title('c^{T}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.8:0.4:1.2);
axis(AX1, [a(1) a(end) 0.8 1.2]);
set(P,NameArray,ValueArray);
subplot(5,2,7);P = plot(a,squeeze(w3(:,lowb)),a,squeeze(w3(:,hghb)));
title('wages','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',1:1:4);
axis(AX1, [a(1) a(end) 1 4]);
set(P,NameArray,ValueArray);
%version 1: with bond spreads
subplot(5,2,8);P = plot(a,squeeze(sp3(:,lowb)),a,squeeze(sp3(:,hghb)));
title('spreads','Fontweight','normal');
set(P,NameArray0,ValueArray0);ylabel('%');
AX1 = gca;box on;
set(AX1,'YTick',0:5:10);
axis(AX1, [a(1) a(end) -0.4 10]);
%version 2: with current account
subplot(5,2,9);P = plot(a,-squeeze(ca3(:,lowb)),a,-squeeze(ca3(:,hghb)));
%xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
title('CA balance','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-0.1:0.1:0.3);
axis(AX1, [a(1) a(end) -0.1 0.3]);
set(P,NameArray0,ValueArray0);
%version 3: with (net) taxes
subplot(5,2,10);P = plot(a,squeeze(tau3(:,lowb)),a,squeeze(tau3(:,hghb)));
set(P,NameArray,ValueArray);
hold on;
PP = plot(a,squeeze(nettax3(:,lowb)),a,squeeze(nettax3(:,hghb)));
title('(net) taxes','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.2:0.5:1.2);
axis(AX1, [a(1) a(end) 0.2 1.2]);
set(PP,NameArray,ValueArray5);
%set(P,NameArray0,ValueArray0);
%subplot(6,2,11);P = plot(b,squeeze(fiscal2(lowa,1,:)),b,squeeze(fiscala2(lowa,1,:)),b,squeeze(fiscalr2(hgha,1,:)),b,squeeze(fiscala2(hgha,1,:)));
%title('primary surplus'); %xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
%subplot(4,2,8);P = plot(b,squeeze(taur2(lowa,1,:)),b,squeeze(taua2(lowa,1,:)),b,squeeze(taur2(hgha,1,:)),b,squeeze(taua2(hgha,1,:)));
%title('T'); xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
%set(P,NameArray,ValueArray);



%% Plot variables as function of debt for two yT levels
hr2  = h2;
cnr2 = cn2;
gnr2 = gn2;
ctr2 = ct2;
pr2  = p2;
taur2= tau2;
nettaxr2= nettax2;
fiscalr2  = fiscal2;
nfiscalr2 = nfiscal2;
bpr2 = bp2;
wwr2 = ww2;
qr2  = q2;
spr2 = 100*((1+qr2.^(-1)-lambda)./(1+r)-1);
car2 = ca2;

%optimal policies and prices: repayment vs autarky
%NameArray = {'LineStyle','Color','LineWidth'};
%ValueArray = {'-',':','-',':';'red','red','blue','blue';1.5,1.5,1.5,1.5}';

%NameArray0 = {'LineStyle','Color','LineWidth'};
%ValueArray0 = {'-','-';'red','blue';1.5,1.5}';
NameArray = {'LineStyle','Color','LineWidth'};
ValueArray = {'--','--','-','-';[1  0.4  0.4],[1  0.4  0.4],[0.078  0.1686  0.549],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray2 = {'--','-','--','-';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray4 = {'--','--','-','-';'red','blue','red','blue';1.5,1.5,1.5,1.5}';
%ValueArray5 = {':',':',':',':';'red','blue','red','blue';1.5,1.5,1.5,1.5}';
%ValueArray5 = {':',':',':',':';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray5 = {':',':',':',':';[1  0.4  0.4],[1  0.4  0.4],[0.078  0.1686  0.549],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';

NameArray0 = {'LineStyle','Color','LineWidth'};
ValueArray0 = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';

[ilowa_aut]=find(def2(lowa,1,:)==0);
[ilowa_rep]=find(def2(lowa,1,:)==1);
hr2(lowa,1,ilowa_aut) = NaN;
ha2(lowa,1,1:nb) = haut2(lowa,1);
ha2(lowa,1,ilowa_rep) = NaN;

[ihgha_aut]=find(def2(hgha,1,:)==0);
[ihgha_rep]=find(def2(hgha,1,:)==1);
hr2(hgha,1,ihgha_aut) = NaN;
ha2(hgha,1,1:nb) = haut2(hgha,1);
ha2(hgha,1,ihgha_rep) = NaN;

pr2(lowa,1,ilowa_aut) = NaN;
pa2(lowa,1,1:nb) = paut2(lowa,1);
pa2(lowa,1,ilowa_rep) = NaN;

pr2(hgha,1,ihgha_aut) = NaN;
pa2(hgha,1,1:nb) = paut2(hgha,1);
pa2(hgha,1,ihgha_rep) = NaN;

gnr2(lowa,1,ilowa_aut) = NaN;
gna2(lowa,1,1:nb) = gnaut2(lowa,1);
gna2(lowa,1,ilowa_rep) = NaN;

gnr2(hgha,1,ihgha_aut) = NaN;
gna2(hgha,1,1:nb) = gnaut2(hgha,1);
gna2(hgha,1,ihgha_rep) = NaN;

bpr2(lowa,1,ilowa_aut) = NaN;
bpr2(hgha,1,ihgha_aut) = NaN;

cnr2(lowa,1,ilowa_aut) = NaN;
cna2(lowa,1,1:nb) = cnaut2(lowa,1);
cna2(lowa,1,ilowa_rep) = NaN;

cnr2(hgha,1,ihgha_aut) = NaN;
cna2(hgha,1,1:nb) = cnaut2(hgha,1);
cna2(hgha,1,ihgha_rep) = NaN;

ctr2(lowa,1,ilowa_aut) = NaN;
cta2(lowa,1,1:nb) = ctaut2(lowa,1);
cta2(lowa,1,ilowa_rep) = NaN;

ctr2(hgha,1,ihgha_aut) = NaN;
cta2(hgha,1,1:nb) = ctaut2(hgha,1);
cta2(hgha,1,ihgha_rep) = NaN;

wwr2(lowa,1,ilowa_aut) = NaN;
wwa2(lowa,1,1:nb) = wwaut2(lowa,1);
wwa2(lowa,1,ilowa_rep) = NaN;

wwr2(hgha,1,ihgha_aut) = NaN;
wwa2(hgha,1,1:nb) = wwaut2(hgha,1);
wwa2(hgha,1,ihgha_rep) = NaN;

taur2(lowa,1,ilowa_aut) = NaN;
taua2(lowa,1,1:nb) = tauaut2(lowa,1);
taua2(lowa,1,ilowa_rep) = NaN;

taur2(hgha,1,ihgha_aut) = NaN;
taua2(hgha,1,1:nb) = tauaut2(hgha,1);
taua2(hgha,1,ihgha_rep) = NaN;

nettaxr2(lowa,1,ilowa_aut) = NaN;
nettaxa2(lowa,1,1:nb) = nettaxaut2(lowa,1);
nettaxa2(lowa,1,ilowa_rep) = NaN;

nettaxr2(hgha,1,ihgha_aut) = NaN;
nettaxa2(hgha,1,1:nb) = nettaxaut2(hgha,1);
nettaxa2(hgha,1,ihgha_rep) = NaN;

fiscalr2(lowa,1,ilowa_aut) = NaN;
fiscala2(lowa,1,1:nb) = fiscalaut2(lowa,1);
fiscala2(lowa,1,ilowa_rep) = NaN;

fiscalr2(hgha,1,ihgha_aut) = NaN;
fiscala2(hgha,1,1:nb) = fiscalaut2(hgha,1);
fiscala2(hgha,1,ihgha_rep) = NaN;

nfiscalr2(lowa,1,ilowa_aut) = NaN;
nfiscala2(lowa,1,1:nb) = nfiscalaut2(lowa,1);
nfiscala2(lowa,1,ilowa_rep) = NaN;

nfiscalr2(hgha,1,ihgha_aut) = NaN;
nfiscala2(hgha,1,1:nb) = nfiscalaut2(hgha,1);
nfiscala2(hgha,1,ihgha_rep) = NaN;

car2(lowa,1,ilowa_aut) = NaN;
car2(hgha,1,ihgha_aut) = NaN;

qr2(lowa,1,ilowa_aut) = NaN;
qr2(hgha,1,ihgha_aut) = NaN;

spr2(lowa,1,ilowa_aut) = NaN;
spr2(hgha,1,ihgha_aut) = NaN;

endb = 2;

figure;
subplot(5,2,1);P = plot(b,100*squeeze(1-hr2(lowa,1,:)),b,100*squeeze(1-ha2(lowa,1,:)),b,100*squeeze(1-hr2(hgha,1,:)),b,100*squeeze(1-ha2(hgha,1,:)));
title('unemployment','Fontweight','normal');
AX1 = gca;box on;ylabel('%');
set(AX1,'YTick',0:10:20);
axis(AX1, [min(b) endb -1 20]);
set(P,NameArray,ValueArray);
subplot(5,2,2);P = plot(b,squeeze(pr2(lowa,1,:)),b,squeeze(pr2(hgha,1,:)),b,squeeze(pa2(lowa,1,:)),b,squeeze(pa2(hgha,1,:)));
title('p^{N}','Fontweight','normal'); 
%lgd=legend('low y^T','high y^T');lgd.Location='SouthEast';
legend('low y^T','high y^T','Location','SouthWest','Orientation','vertical');
AX1 = gca;box on;
set(AX1,'YTick',1:3:7);
axis(AX1, [min(b) endb 1 7]);
set(P,NameArray,ValueArray2);
subplot(5,2,3);P = plot(b,squeeze(gnr2(lowa,1,:)),b,squeeze(gna2(lowa,1,:)),b,squeeze(gnr2(hgha,1,:)),b,squeeze(gna2(hgha,1,:)));
title('g^{N}','Fontweight','normal'); 
AX1 = gca;box on;
set(AX1,'YTick',0:0.2:0.4);
axis(AX1, [min(b) endb 0 0.4]);
set(P,NameArray,ValueArray);
%subplot(5,2,4);P = plot(b,squeeze(bpr2(lowa,1,:)),b,squeeze(bpr2(hgha,1,:)));
subplot(5,2,4);P = plot(b,-squeeze(bpr2(lowa,1,:))/(lambda+r),b,-squeeze(bpr2(hgha,1,:))/(lambda+r));
title('debt','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-3:3:3);
axis(AX1, [min(b) endb -3 3]);
set(P,NameArray0,ValueArray0);
subplot(5,2,5);P = plot(b,squeeze(cnr2(lowa,1,:)),b,squeeze(cna2(lowa,1,:)),b,squeeze(cnr2(hgha,1,:)),b,squeeze(cna2(hgha,1,:)));
title('c^{N}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.8:0.1:1);
axis(AX1, [min(b) endb 0.8 1]);
set(P,NameArray,ValueArray);
subplot(5,2,6);P = plot(b,squeeze(ctr2(lowa,1,:)),b,squeeze(cta2(lowa,1,:)),b,squeeze(ctr2(hgha,1,:)),b,squeeze(cta2(hgha,1,:)));
title('c^{T}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.5:0.5:1.5);
axis(AX1, [min(b) endb 0.5 1.5]);
set(P,NameArray,ValueArray);
subplot(5,2,7);P = plot(b,squeeze(wwr2(lowa,1,:)),b,squeeze(wwa2(lowa,1,:)),b,squeeze(wwr2(hgha,1,:)),b,squeeze(wwa2(hgha,1,:)));
title('wages','Fontweight','normal');xlim([min(b) max(b)]);
AX1 = gca;box on;
set(AX1,'YTick',0:3:6);
axis(AX1, [min(b) endb 0 6]);
set(P,NameArray,ValueArray);
%version 1: with bond spreads
subplot(5,2,8);P = plot(b,squeeze(spr2(lowa,1,:)),b,squeeze(spr2(hgha,1,:)));
title('spreads','Fontweight','normal');
set(P,NameArray0,ValueArray0);ylabel('%');
AX1 = gca;box on;
set(AX1,'YTick',0:25:50);
axis(AX1, [min(b) endb -2 50]);
%version 2: with current account
subplot(5,2,9);P = plot(b,-squeeze(car2(lowa,1,:)),b,-squeeze(car2(hgha,1,:)));
%xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
title('CA balance','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-0.5:0.5:0.5);
axis(AX1, [min(b) endb -0.5 0.5]);
set(P,NameArray0,ValueArray0);
%version 3: with (net) taxes
subplot(5,2,10);P = plot(b,squeeze(taur2(lowa,1,:)),b,squeeze(taua2(lowa,1,:)),b,squeeze(taur2(hgha,1,:)),b,squeeze(taua2(hgha,1,:)));
set(P,NameArray,ValueArray);
hold on;
PP = plot(b,squeeze(nettaxr2(lowa,1,:)),b,squeeze(nettaxa2(lowa,1,:)),b,squeeze(nettaxr2(hgha,1,:)),b,squeeze(nettaxa2(hgha,1,:)));
%xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
title('(net) taxes','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0:0.5:1.5);
axis(AX1, [min(b) endb 0 1.5]);
set(PP,NameArray,ValueArray5);
%set(P,NameArray0,ValueArray0);
%subplot(6,2,11);P = plot(b,squeeze(fiscal2(lowa,1,:)),b,squeeze(fiscala2(lowa,1,:)),b,squeeze(fiscalr2(hgha,1,:)),b,squeeze(fiscala2(hgha,1,:)));
%title('primary surplus'); %xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
%subplot(4,2,8);P = plot(b,squeeze(taur2(lowa,1,:)),b,squeeze(taua2(lowa,1,:)),b,squeeze(taur2(hgha,1,:)),b,squeeze(taua2(hgha,1,:)));
%title('T'); xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
%set(P,NameArray,ValueArray);

dbstop;


%% plot time series
seriesp   = load('seriesp.txt');
seriesr   = load('seriesr.txt');
seriesh   = load('seriesh.txt');
seriesb   = load('seriesb.txt');
seriesct  = load('seriesct.txt');
seriesm   = 1; %load('seriesm.txt');
seriesgn  = load('seriesgn.txt');
seriestau = load('seriestau.txt');    
seriesa   = load('seriesa.txt');    
seriesbor = load('seriesbor.txt');    

n = length(seriesp);

seriescn= seriesh.^alpha-seriesgn;
seriesy = seriesa.*seriesm.^gamma+seriesp.*seriesh.^alpha; % -pm*seriesm;
seriesgny = seriesgn./seriesy;

time =1:n;

ValueArray1 = {'-';'blue';1.5}';
ValueArray2 = {'-','--';'red','blue';1.5,1.5}';
ValueArray3 = {'-','--','-';'red','blue','black';1.5,1.5,1.5}';

% No shaded areas for default episodes
% figure;
% subplot(3,2,1); P=plot(time,seriesa');title('a');
% set(P,NameArray,ValueArray1);
% subplot(3,2,2);P=plot(time,seriesp);title('p');
% set(P,NameArray,ValueArray1);
% subplot(3,2,3);P=plot(time,seriesct);
% %AX=legend('c_{T}','Location','NorthWest','Orientation','Horizontal');
% title('c_T');
% set(P,NameArray,ValueArray1);
% %LEG = findobj(AX,'type','text');
% %set(LEG,'FontSize',7);
% subplot(3,2,4); PPP=plot(time,seriescn,time,seriesh,time,seriesgn);
% axis([1 n -0.05 1.05]);
% AX=legend('h','c_{N}','g_{N}','Location','NorthWest','Orientation','Horizontal');title('h,c_{N} & g_{N}');
% set(PPP,NameArray,ValueArray3);
% LEG = findobj(AX,'type','text');
% set(LEG,'FontSize',7);
% subplot(3,2,5); P=plot(time,seriesr);title('spreads');
% set(P,NameArray,ValueArray1);axis([1 n -0.01 0.1]);
% subplot(3,2,6); P=plot(time,-seriesb);title('b');
% set(P,NameArray,ValueArray1);


seriesbor = 1-seriesbor;
seriesbor(199)=1;
seriesbor(200)=0;
seriesh(200) = 0;
seriesr(200) = 0.1;
ntime = 200;

figure;
subplot(3,2,1); 
shadedTimeSeries(time,seriesa',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('a');
AX1 = gca;box on;
set(AX1,'YTick',0.85:0.1:1.15);
axis(AX1, [time(1) time(ntime) 0.85 1.15]);
subplot(3,2,2);
shadedTimeSeries(time,seriesp',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('p');
AX1 = gca;box on;
set(AX1,'YTick',0:0.2:1);
axis(AX1, [time(1) time(ntime) 0.95*min(seriesp) 1.05*max(seriesp)]);
subplot(3,2,3);
shadedTimeSeries(time,seriesct',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('c_T');
AX1 = gca;box on;
set(AX1,'YTick',0.05:0.05:1.15);
axis(AX1, [time(1) time(ntime) 0.9*min(seriesct)  1.1*max(seriesct)]); 
subplot(3,2,4); 
shadedTimeSeries(time,seriesh',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('h,c_N & g_N');
AX1 = gca;box on;
axis(AX1, [time(1) time(ntime) 0 1.05]);
hold on;
H0 = line(time,seriesh','Color','b','LineWidth',2);
H1 = line(time,seriescn','Color','r','LineWidth',2,'LineStyle','--');
H2 = line(time,seriesgn','Color','g','LineWidth',2,'LineStyle',':');
hh = [H0, H1,H2];
legend(hh,'h','c_N','g_N','Location','SouthEast','Orientation','Horizontal');
AX1 = gca;box on;
axis(AX1, [time(1) time(ntime) 0 1.05]);
subplot(3,2,5); 
shadedTimeSeries(time,seriesr',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('spreads');
AX1 = gca;box on;
set(AX1,'YTick',0.0:0.025:0.1);
axis(AX1, [time(1) time(ntime) 0 0.1]);
subplot(3,2,6); P=plot(time,-seriesb);title('b');
shadedTimeSeries(time,-seriesb',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('b');
AX1 = gca;box on;
%set(AX1,'YTick',0.0:0.02:ceil(max(-1.1*seriesb/0.001))*0.001);
axis(AX1, [time(1) time(ntime) 0 max(-seriesb)]);

dbstop;